First-Principles Description of Correlation Effects in Layered Materials 
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We present a first-principles description of anisotropic materials characterized by having both 
weak (dispersion-like) and strong covalent bonds, based on the Adiabatic-Connection Fluctuation- 
Dissipation Theorem within Density Functional Theory. For hexagonal boron nitride the in-plane 
and out of plane bonding as well as vibrational dynamics are well described both at equilibrium 
and when the layers are pulled apart. Also bonding in covalent and ionic solids is described. The 
formalism allows to ping-down the deficiencies of common exchange-correlation functionals and 
provides insight towards the inclusion of dispersion interactions into the correlation functional. 

PACS numbers: 71.10.-w, 71.15.Mb, 71.15.Nc, 63.20.-e 
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Soft matter, biomolecules, and layered materials are 
examples of weakly bonded compounds, which consti- 
tute stringent tests for ab-initio electronic structure cal- 
culations. A successful description of such systems re- 
quires an accurate treatment of a whole range of inter- 
actions, from short-range covalent to long-range van der 
Waals (vdW) forces. Moreover, their study provides valu- 
able hints about how to devise new approximations that 
describe at a reasonable computational cost very differ- 
ent bonding regimes. In this realm, layered materials 
represent a unique case. The Kohn-Sham (KS) method 
under the simple local density approximation (LDA) [l| 
is known to give reasonable structural properties near the 
equilibrium interlayer distance of systems like graphite Q 
and hexagonal boron nitride (/i-BN) 0,0- Surprisingly 
the more elaborated exchange-correlation (XC) energy 
functional based on the generalized-gradient approxima- 
tion (GGA) [3, while provides a good in-plane bonding 
fails dramatically in the description of the equilibrium 
interlayer distance in /i-BN and graphite, leading to com- 
pletely unphysical results 0, 0] ■ The mechanism respon- 
sible for this poor performance is not specific to GGA 
but to all functionals based on exact-exchange plus 
some LDA correlation (EXXc) (as illustrated below). 

It has been hypothesized that this failure is due to the 
lack of description of vdW forces, which are the manifes- 
tation of long-range correlation effects fTJ . However LDA, 
without accounting for vdW interactions, does bind the 
/i-BN and graphite layers. More importantly, recently 
proposed approximate non-local prescriptions beyond 
LDA or GGA in the vdW Density Functional (vdW-DF) 
framework |||, while describing vdW forces, result in a 
large overestimation of the equilibrium layer-layer dis- 
tance Q. Thus, available XC functionals cannot de- 
scribe correctly a layered material at different geomet- 



rical configurations. Consequently, there is a basic need 
to understand this severe breakdown and how to define 
a stable, successful scheme to calculate such structural 
properties. This is the goal of the present Letter. We 
demonstrate that a fully microscopic, parameter-free ap- 
proach based on the Adiabatic-Connection Fluctuation- 
Dissipation Theorem (ACFDT) [10j correctly describes 
the mechanical properties of /i-BN near and far from the 
equilibrium configuration. We show that ACDFT does 
bind the /i-BN layers at an equilibrium distance in excel- 
lent agreement with the experiment, accounting by con- 
struction for vdW forces In contrast, while LDA 
overbinds /i-BN, GGA, EXXc and, obviously, any other 
prescription with a local or semi-local correlation part 
fail in describing the layers bonding. 



/i-BN has a wealth of important applications for high 
resistance and inert nanodevices |J 0, 0] . Further- 
more, the ft-BN layer-layer interaction is weaker than 
in graphite (the paradigmatic example of layered ma- 
terial), thus being a better scenario to assess and un- 
derstand the performance of KS-based schemes. The 
ACFDT has been a ppli ed so far only to jellium sys- 
tems , molecules |l4( , or compact solids . How- 
ever it starts to be consider a promising ab-initio method 
because: i) it provides a reliable first-principles descrip- 
tion of XC effects at a lower numerical cost than Quan- 
tum Chemistry methods or statistical Quantum Monte- 
Carlo techniques; ii) it is naturally built into the KS 
scheme as a way to construct advanced XC energy func- 
tionals depending on the electron density n (r) and on 
the set of occupied [4> v (r)] and empty [<f> c (r)] KS or- 
bitals. Indeed, according to the ACFDT, the total XC 
energy comprises of the exact-exchange energy Ex — 
-3 / drdr' |£° cc 0* (r) fa (r')| 2 / |r - r'| , plus the corre- 
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lation energy given exactly by 
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dX Tr (w [x\ (iu) - Xo (iu)}) , (1) 



(we will use atomic units unless otherwise stated). Tr de- 
notes the spacial trace, and x\ (iu) is the density-density 
response function at imaginary frequencies of a fictitious 
electron-system interacting through a scaled Coulomb 
potential Aw (r) = X/r whose ground-state density equals 
the actual one. Thus, xo (iu) is the non-interacting re- 
sponse function of this fictitious KS system: 



Xo (r,r';iu) = 25ft 



<j>% ( r) <j> v (rQ# (r')</> e (r) 
iu + (e v - e c ) 



(2) 



where e v (e c ) are the occupied (empty) KS eigenener- 
gies. The interacting response x can be obtained in the 
framework of time-dependent Density Functional The- 
ory (TDDFT) HI, [13 by solving a Dyson-like equation: 

X\ (iu) = (l - Xo (iu) Xw + /xaA (iu) ) x\ (iu) , (3) 

where /xc,a (iu) is the XC kernel of the fictitious system 
with the scaled interaction Xw. The Random Phase Ap- 
proximation (RPA) to x is obtained by setting fxc = 0. 
Note that, in contrast to other procedures put forward 
in the literature 0, @ , this scheme provides a correlation 
functional built at the same level than exact exchange. 

The XC effects introduced by fxc on x for the the cal- 
culation of Ec are considerably weaker than for the de- 
scription of the optical properties of /i-BN |l8|, as y(iu) 
is a very smooth function on the imaginary axis |l9| . 
Therefore, we set /xc,a (iu) = in Eq. J3J) and correct 
the well known deficiency of the RPA total correlation 
energy via a LDA-like correction term (a prescrip- 
tion often labeled as RPA+) |2l|. This scheme will be 
denoted in what follows as EXX/RPA+. Moreover, in- 
stead of solving selfconsistently the KS equations for the 
EXX/RPA+ potential H3, we evaluated the XC energy 
using the KS-LDA electronic structure as commonly done 
in Ab-initio Many-Body calculations 01- We have ex- 
plicitly checked that the use of GGA or EXXc KS states 
as a basis for the EXX/RPA+ expressions has little effect 
(less than 2 %) on the final XC energy. 

Some technical details are relevant. The KS response 
function is computed in reciprocal space using well con- 
verged KS orbitals and eigenenergies obtained from the 
ABINIT code |23. The interacting responses xa are then 
obtained for each vector in the Brillouin zone (BZ), imag- 
inary frequency u, and coupling parameter A using the 
self code |24|. Finally, the RPA correlation energy 
is calculated according to Eq. @J using a sixth order 
Gauss-Legendre (GL) sampling for A integration. The 
imaginary frequency integrals are done using two GL 
grids, each typically comprising 24 points. This allows 
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FIG. 1: (Color on-line) The LDA, EXX/RPA+, and GGA 
exchange and correlation energies per unit cell for A-BN as a 
function of the layer-layer distance d (the correlation energy 
of Si as a function of the lattice constant a is showed in the 
inset). 



a well converged evaluation of the contributions from 
both small and large imaginary frequencies. The con- 
vergence with respect to the BZ sampling has been also 
carefully checked (the exchange energy is very sensitive 
to this sampling). Nonetheless, the critical convergence 
parameters from a computational point of view are the 
number of unoccupied bands and reciprocal space vec- 
tors in Eq. (J2J . The use of a plane- wave representation 
allows a systematic simultaneous convergence with re- 
spect to these two parameters. Further technical details, 
together with more applications to metallic systems and 
graphite, will be presented elsewhere. 

Before studying the weak interlayer interaction in 
/i-BN, it is illustrative to analyze the impact of the 
exchange-correlation effects at the EXX/RPA+ level for 
two prototype compact structures: bulk Si (as covalent 
semiconductor) and NaCl (as ionic material). It is known 
that the standard XC functionals, like GGA, provide a 
good description of those systems. This is confirmed in 
the present ACFDT scheme as the impact in the lat- 
tice constant of the EXX/RPA+ XC energy correction 
is below 0.5%. However, the hybrid EXXc overcorrects 
the GGA results specially in the low density regime of 
bulk Si due to the LDA overestimation of the correla- 
tion energy for low densities (see inset of FigQ]). Nev- 
ertheless, this does not compromise very much the ac- 
curacy of the EXXc approach since in this case the ex- 
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FIG. 2: (Color on-line) Variation of the total energy per 
unit cell of ft-BN versus the layer-layer distance d. While 
EXX/RPA+ binds the ft-BN layers, GGA and EXXc dra- 
matically fail. EXX/PRA+ improves the equilibrium inter- 
layer distance compared to vdW-DF and LDA, whose good 
performance is merely due to the cancellation of errors. The 
EXX/RPA+ error bars come from the extrapolation tech- 
nique used to extract the correlation energy from values cor- 
responding to finite number of bands and reciprocal space 
calculations. The experimental equilibrium distance is taken 
from Ref. 



change energy dominates (Ex (11) —-Ex (9-6) s=y 0.2 Ha ^> 
E c (11) - E c (9.6) w 0.01 Ha) jig. 

The previous scenario, where the variations of the ex- 
change energy is at least one order of magnitude larger 
than correlation, changes completely in the case of ft- 
BN. As displayed in Fig.^ the variations of Ex and Eq 
around the equilibrium interlayer distance d eq in ft-BN 
are of the same order of magnitude (few mHa) and at 
least one order of magnitude smaller than for Si or NaCl. 
Thus, small errors are more clearly visible and correlation 
plays a much more important role than in the case of the 
compact structures. Hence, a consistent description of 
exchange and correlation is needed to reproduce the deli- 
cate energetic interplay in the weak interlayer bonding of 
ft-BN. This consistency is absent in the EXXc prescrip- 
tion and, as a consequence, it provides too large (or even 
unbound) interlayer equilibrium distances, as showed in 
Fig-El To a lesser extent, this problem is also present 
in the vdW-DF scheme designed to treat vdW forces in 
an approximate manner that largely overestimates the 
equilibrium distance (6.86 compared to 6.24 a.u.), while 
underestimating the bulk modulus. The GGA follows 
the very same trend than the EXXc. This puts in evi- 
dence the fact that GGA exchange greatly improves upon 
the LDA (and, hence, provides better chemical properties 
of molecules), but does not do a so good job for corre- 
lation. Instead EXX/RPA+, which treats consistently 
exchange and correlation, correctly predicts ft-BN to be 



stable, with structural parameters in excellent agreement 
with experiment (see Fig.|2J). Furthermore, in contrast to 
LDA, GGA, and EXXc, EXX/RPA+ naturally embod- 
ies the correct long-range d~ vdW power law depen- 
dence of the energy Q as a function of interlayer dis- 
tance, thus merging continuously the covalent and vdW 
regimes. This is confirmed in Fig. [2] where, for d larger 
than 7.5 a.u., we observe that LDA deviates from the cor- 
rect total energy dependence (given by EXX/RPA+). By 
looking at the LDA XC contributions in Fig. ^ we see 
that even if both Ex and Ec are poorly described, the 
XC energy dependence on the interlayer distance is sim- 
ilar to the ACDFT result. Nevertheless, the slope of the 
LDA XC energy is larger than the ACDFT one. This dif- 
ference is responsible for the smaller interplane distance 
in the LDA as compared to EXX/RPA+, being the latter 
(6.25 a.u.) in stringent ag reement with the experimental 
value c£P = 6.24 a.u. Eg. Moreover in the ACFDT the 
LDA bulk modulus B is increased by about 5% in better 
agreement with the experimental value (in contrast to 
the vdW-DF scheme [7J where B is reduced). Therefore, 
while the relative success of the LDA around the equi- 
librium configuration is merely due to a cancellation of 
errors the EXX/RPA+ scheme provides better structural 
properties describing equally well the correlation and the 
exchange energy. This shows that the ACFDT Ab-initio 
scheme is suitable to study very different regimes, like Si 
and NaCl, where the exchange energy dominates, or like 
ft-BN, where the success is bound to a correct balance of 
correlation and exchange contributions. 

An additional severe test of the accuracy of the 
EXX/RPA+ method is given by the calculation of the 
vibrational spectra of ft-BN. LDA and GGA in-plane 
lattice dynamics Q reproduces fairly well the available 
experimental results. Similarly, the results of a frozen- 
phonon calculation within the present EXX/RPA+ 
scheme give in-plane modes at the high-symmetry points 
within few percents from experiments (and LDA or 
GGA) . The situation changes considerably for interplane 
modes where the LDA and GGA phonons differ by up to 
30 % when calculated at the experimental interplane dis- 
tance (this is the only way to compare the results as GGA 
does not bind ft-BN and LDA overbinds). Indeed, modes 
related to interlayer vibrations would be very much sensi- 
tive to errors in the correlation functional (as the typical 
variations in the correlation energy that has to be re- 
solved are up to two orders of magnitude smaller for inter- 
layer phonons than for structural characterization). Fur- 
thermore, recent second-order Raman experiments |26j | 
have suggested that the interplane modes are greatly un- 
derestimated in the theoretical calculations. In partic- 
ular they quote a frequency value for the phonon 
mode (where two BN planes rigidly slide against each 
other) of 310 cm -1 , which is about a factor of three larger 
than the LDA (90 cm -1 ) or GGA (122 cm" 1 ) values^, 
computed at the experimental lattice parameters. If the 
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experimental analysis based on overtones was correct it 
would imply that correlation effects (dispersion forces) 
are very important for those vibrations. This point can 
be directly addressed by the present EXX/RPA+ scheme 
as we have shown that it describes properly the struc- 
tural properties of ft-BN and accounts at the same time 
by construction for the long-range vdW forces. Surpris- 
ingly the calculated EXX/RPA+ phonon mode en- 
ergy is 130±10 cm -1 , regardless the use of the theoretical 
or experimental lattice parameters. Therefore, it is clear 
that improving the description of XC effects upon LDA 
brings the T± stiffer, still not close to the experimental 
assignment. However, this huge discrepancy can be dras- 
tically reduced if in the analysis of the experimental data 
we allow for the whole set of double-phonon processes 
(sum and difference) consistent with the selection rules 
of Raman scattering instead of considering only phonon- 
overtones. If we do so, the new frequency assignment 
for the r^" phonon mode is about 125 cm -1 |2jj, still 
in disagreement with LDA value but very close to the 
EXX/RPA+ result and to the GGA value (in spite of the 
fact that it does not bind h-BN). Therefore, EXX/RPA+ 
solves unambiguously this puzzling situation providing 
an accurate method to evaluate the role played by the 
exchange-correlation effects. 

In summary, we have shown that EXX/RPA+ cor- 
rectly describes the structural properties, the on-plane 
and the out-of-plane phonon modes of /i-BN. While both 
GGA and EXXc do not predict the ft-BN layers bond- 
ing and LDA underestimates both the interlayer phonon 
mode frequency and interlayer distance, EXX/RPA+ 
provides a consistent, accurate and successful alterna- 
tive scheme. Thus, we have seed light about the per- 
formance of standard ab-initio techniques when applied 
to layered systems where any improvements upon LDA 
require equally accurate treatments of exchange and cor- 
relation effects. The EXX/RPA+ scheme presented here 
opens the field to more elaborated ab-initio applications 
in the wide field of sparse systems and soft-matter. 
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